In [264]:

    

using PyPlot
using PyCall
@pyimport matplotlib.patches as patches

function block(A,s1=5,s2=5)  
    m, n=size(A)
    mi=[i:min(i+s1-1,m) for i=1:s1:m]
    ni=[j:min(j+s2-1,n) for j=1:s2:n]
    [A[i,j] for i in mi, j in ni]
end

unblock(A)=hvcat(tuple([size(A,2) for i=1:size(A,1)]...), A'...) 

# Key BLAS3 block operator
# A few block apply Q's
import Base.LinAlg.QRPackedQ, Base.LinAlg.Ac_mul_B, Base.*
Ac_mul_B(Q::Base.LinAlg.QRPackedQ,B)=hcat([block(Q'*unblock(B[:,j])) for j in 1:size(B,2)]...)        
*{T<:Matrix}(A::Array{T},Q::Base.LinAlg.QRPackedQ)=vcat([block(unblock(A[i,:])*Q) for i in 1:size(A,1)]...)   
                                        
# A few block qrQ's
qrQ{T<:Number}(A::Array{T})=qrfact(A)[:Q]
qrQ{T<:Matrix}(A::Array{T})=qrQ(unblock(A))  
qrQc{T<:Matrix}(A::Array{T})=qrQ(unblock(A)')
                                 
spyt{T<:Number}(A::Array{T})=spy(A, precision=sqrt(eps(T)), alpha=0.8)                                    
spy{T<:Matrix}(A::Array{T})=spy(unblock(A), precision=sqrt(eps()), alpha=0.7)  # Threshold for convenience

#Convenience drawing functions
rect_red(x0, y0, dx, dy)=gca()[:add_patch](patches.Rectangle((x0-1.5,y0-1.5),dx,dy,linewidth=5,facecolor="none",edgecolor="red"))
rect_blue(x0, y0, dx, dy)=gca()[:add_patch](patches.Rectangle((x0-1.5,y0-1.5),dx,dy,linewidth=3,facecolor="none",edgecolor="blue"))
rect_gray(x0, y0, dx, dy)=gca()[:add_patch](patches.Rectangle((x0-1.5,y0-1.5),dx,dy,linewidth=1,facecolor="none",edgecolor="#aaaaaa"))
rect_red(yr, xr)=rect_red(xr.start,yr.start,xr.len,yr.len)
rect_blue(yr, xr)=rect_blue(xr.start,yr.start,xr.len,yr.len)
rect_gray(yr, xr)=rect_gray(xr.start,yr.start,xr.len,yr.len)


    

    
Out[264]:





rect_gray (generic function with 2 methods)

In [265]:

    

N=15
s=5; #block size
A=block(randn(N,N));
Q=qrQ(A[1,1])
A[1,:]=Q'*A[1,:]

spy(A)
rect_red(1:s,1:s)  #QR from this block
rect_blue(1:s,1:N) #QR applied to this block


    

    

















    
Out[265]:





PyObject <matplotlib.patches.Rectangle object at 0x4959a210>

In [266]:

    

for j=2:3
    Q=qrQ(A[[1,j],1])        # There is a special QR for a "domino"=[◥;◼]
    A[[1,j],:]=Q'*A[[1,j],:] # Q' acts like a square matrix even though Q was "thin" 
    rect_red(1,(j-1)*s+1,s,s)#QR from this block
end 
spy(A)


    

    

















    
Out[266]:





PyObject <matplotlib.image.AxesImage object at 0x5988ebd0>

In [267]:

    

#Now do this from the right
Q=qrfact(A[1,2]')[:Q]  # We really need an LQ
A[:,2]*=Q
spy(A)
rect_red(s+1,1,s,s)  #QR from this block
rect_blue(1:s,s+1:N) #QR applied to this block


    

    

















    
Out[267]:





PyObject <matplotlib.patches.Rectangle object at 0x541aecd0>

In [268]:

    

for j=3:3
    Q=qrQ(unblock(A[1,[2,j]])') # This is a special LQ for a "domino"=[◥;◼]'                       
    A[:,[2,j]]*=Q
    rect_red((j-1)*s+1,1,s,s)#QR from this block
end
spy(A)


    

    

















    
Out[268]:





PyObject <matplotlib.image.AxesImage object at 0x4d1ea150>

In [269]:

    

function BandBidiagonal(A::Array{Array{Float64,2},2}) 
    (m,n)=size(A)
    for i=1:n
        Q=qrQ(A[i,i])
        A[i,:]=Q'*A[i,:]
        for j=(i+1):n
          Q=qrQ(A[[i,j],i])         # QR for a "domino"=[◥;◼]
          A[[i,j],:]=Q'*A[[i,j],:]  # Q' acts like a square matrix even though Q was "thin" 
        end
        i==n && return A
        Q=qrQ(A[i,i+1]')            # We really need an LQ
        A[:,i+1]*=Q
        for j=(i+2):n
          Q=qrQc(A[i,[i+1,j]])      # LQ for a "domino"=[◥;◼]'                       
          A[:,[i+1,j]]*=Q                       
        end
    end
end


    

    
Out[269]:





BandBidiagonal (generic function with 1 method)

In [270]:

    

A=block(randn(15,15))
A_after=BandBidiagonal(A)
spy(A_after)
[svdvals(unblock(A))[1:5] svdvals(unblock(A_after))[1:5]]


    

    

















    
Out[270]:





5x2 Array{Float64,2}:
 6.94917  6.94917
 6.10788  6.10788
 5.63312  5.63312
 5.15745  5.15745
 4.64259  4.64259

In [271]:

    

A=unblock(A_after)
Q1=qrQ(A[1,2:6]')
A[1:6,2:6]*=Q1
spyt(A)
rect_red(1:1,2:6)  #QR from this block
rect_blue(1:6,2:6) #QR applied to this block


    

    

















    
Out[271]:





PyObject <matplotlib.patches.Rectangle object at 0x4c572c50>

In [272]:

    

Q2=qrQ(A[2:6,2:2])
A[2:6,2:11]=Q2'*A[2:6,2:11]
spyt(A)
rect_red(2:6,2:2)  #QR from this block
rect_blue(2:6,2:11) #QR applied to this block


    

    

















    
Out[272]:





PyObject <matplotlib.patches.Rectangle object at 0x47fa3e90>

In [273]:

    

Q3=qrQ(A[2,7:11]')#Clear the fill-in from the previous block
A[2:11,7:11]*=Q3
spyt(A)
rect_red(2:2,7:11)  #QR from this block
rect_blue(2:11,7:11) #QR applied to this block


    

    

















    
Out[273]:





PyObject <matplotlib.patches.Rectangle object at 0x47fbad50>

In [274]:

    

A=block(randn(25,25),5,5);
A=BandBidiagonal(A);
s=size(A[1,1],1) #Square block size
n=size(A,1)      #Number of blocks
A=unblock(A)
i=1
Q=qrQ(A[i,i+(1:s)]') #odd
A[i+(0:s),i+(1:s)]*=Q
rect_red(i:i,i+(1:s))  #QR from this block
rect_blue(i+(0:s),i+(1:s)) #QR applied to this block
for i=1:s:((n-2)*s)
    Q=qrQ(A[i+(1:s),i+1:i+1]) #even
    A[i+(1:s),i+(1:2s)]=Q*A[i+(1:s),i+(1:2s)]
    rect_red(i+(1:s),i+1:i+1)  #QR from this block
    rect_blue(i+(1:s),i+(1:2s)) #QR applied to this block
    Q=qrQ(A[i+1,i+s+(1:s)]')  #odd
    A[i+(1:2s),i+s+(1:s)]*=Q
    rect_red(i+1:i+1,i+s+(1:s))  #QR from this block
    rect_blue(i+(1:2s),i+s+(1:s)) #QR applied to this block
end
i=(n-2)*s+1
Q=qrQ(A[i+(1:s),i+1:i+1]) #even
A[i+(1:s),i+1:end]=Q*A[i+(1:s),i+1:end]
rect_red(i+(1:s),i+1:i+1)  #QR from this block
rect_blue(i+(1:s),i+1:(n*s)) #QR applied to this block
Q=qrQ(A[i+1,i+s+1:end]')  #odd
A[i+1:end,i+s+1:end]*=Q
rect_red(i+1:i+1,i+s+1:(n*s))  #QR from this block
rect_blue(i+1:(n*s),i+s+1:(n*s)) #QR applied to this block
i=(n-1)*s+1
Q=qrQ(A[i+1:end,i+1:i+1]) #even
A[i+1:end,i+1:end]=Q*A[i+1:end,i+1:end]
rect_red(i+1:(n*s),i+1:i+1)  #QR from this block
rect_blue(i+1:(n*s),i+1:(n*s)) #QR applied to this block
spyt(A)


    

    

















    
Out[274]:





PyObject <matplotlib.image.AxesImage object at 0x52162790>

In [275]:

    

@pyimport matplotlib.animation as anim

function block_bidiagonalize(M,s1=5,s2=5; record_video=false, video_file="blocksvd.mp4")
    if record_video
        video = anim.FFMpegWriter(fps=6, extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])
        video[:setup](gcf(), video_file)
    end
    A=BandBidiagonal(block(M,s1,s2))
    s=size(A[1,1],1) #Square block size
    n=size(A,1)      #Number of blocks
    A=unblock(A)
    for j=1:n*s-2    #bulge chasing: elimination on row j+1
        endb=min(j+s,n*s)
        Q=qrQ(A[j,j+1:endb]')           #xGBCW1
        A[j:endb,j+1:endb]*=Q
        spyt(A)
        rect_red(j:j,j+1:endb)  #QR from this block
        rect_blue(j:endb,j+1:endb) #QR applied to this block
        record_video && video[:grab_frame]()

        lastb=((n-floor(j/s)-1)*s)+j                 #index of last block
        for i=j:s:lastb
            endb, endbp1=min(i+s,n*s), min(i+2s,n*s) #index of end of block and its neighbor
            Q=qrQ(A[i+1:endb,i+1:i+1])  #xGBCW2
            A[i+1:endb,i+1:endbp1]=Q*A[i+1:endb,i+1:endbp1]
            spyt(A)
            rect_red(i+1:endb,i+1:i+1)  #QR from this block
            rect_blue(i+1:endb,i+1:endbp1) #QR applied to this block
            record_video && video[:grab_frame]()
            i==lastb && break
            Q=qrQ(A[i+1,i+s+1:endbp1]') #xGBCW3
            A[i+1:endbp1,i+s+1:endbp1]*=Q
            spyt(A)
            rect_red(i+1:i+1,i+s+1:endbp1)  #QR from this block
            rect_blue(i+1:endbp1,i+s+1:endbp1) #QR applied to this block
            record_video && video[:grab_frame]()
        end
    
        #Gray out stuff
        clf()
        endb=min(j+s,n*s)
        rect_gray(j:j,j+1:endb)
        rect_gray(j:endb,j+1:endb)
        for i=j:s:lastb
            endb, endbp1=min(i+s,n*s), min(i+2s,n*s)
            rect_gray(i+1:endb,i+1:i+1)
            rect_gray(i+1:endb,i+1:endbp1)
            rect_gray(i+1:i+1,i+s+1:endbp1)
            rect_gray(i+1:endbp1,i+s+1:endbp1)
        end
    end
    clf()
    spyt(A)
    if record_video
        video[:grab_frame]()
        video[:finish]()
    end
    A
end


    

    
Out[275]:





block_bidiagonalize (generic function with 3 methods)

In [276]:

    

function block_svdvals(M,s1=5,s2=5)
    A=block_bidiagonalize(M,s1,s2)
    B=Bidiagonal(diag(A), diag(A,1), true)
    svdvals(B)
end

[block_svdvals(M) svdvals(M)]


    

    

















    
Out[276]:





25x2 Array{Float64,2}:
 9.47961   9.47961 
 9.06957   9.06957 
 8.10747   8.10747 
 7.56554   7.56554 
 6.84574   6.84574 
 6.6843    6.6843  
 6.49103   6.49103 
 6.27427   6.27427 
 5.59487   5.59487 
 5.28317   5.28317 
 5.1632    5.1632  
 4.4533    4.4533  
 4.06357   4.06357 
 3.80315   3.80315 
 3.3645    3.3645  
 2.98093   2.98093 
 2.82718   2.82718 
 2.59303   2.59303 
 2.41564   2.41564 
 1.69782   1.69782 
 1.43786   1.43786 
 1.10352   1.10352 
 0.537197  0.537197
 0.470536  0.470536
 0.133043  0.133043

In [277]:

    

M=randn(25,25)
A=block_bidiagonalize(M,5,5)
spyt(A) #imshow(A, cmap="bwr")


    

    

















    
Out[277]:





PyObject <matplotlib.image.AxesImage object at 0x4cb1c250>

In [285]:

    

;ls -l *.mp4
embed_video(filename)=display("text/html", string("""<video autoplay controls><source src="data:video/x-m4v;base64,""",
                            base64(open(readbytes,filename)),"""" type="video/mp4"></video>"""))
embed_video("blocksvd.mp4")


    

    




bash: no job control in this shell






    




-rw-rw-r-- 1 ec2-user ec2-user   172226 Nov  8 23:27 blocksvd.mp4
-rw-rw-r-- 1 ec2-user ec2-user 17802402 Nov  9 00:37 turkey_svd.mp4
-rw-rw-r-- 1 ec2-user ec2-user    52700 Nov  8 21:47 writer_test.mp4






    

In [ ]: